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I. Project Objectives (per proposal) 

1 . Rewrite our existing finite-element code to permit it to be used in global simulations. 

2. Create a forward model of the dynamics of the global lithosphere, permitting predictions 

of how alternative patterns of mantle convection and alternative plate rheologies 
would be expressed as observable geodetic displacements and in situ stresses. 

3. Apply this model in a set of systematic simulation experiments to determine these 

quantities: 

4. First, experiment with parameterized tractions on the bottoms and edges of plates to 

match geodetic plate velocities. Next, optimize tractions on each plate by an iterative 
method. Finally, vary plate rheology to match geodetic intraplate strain rates and in 
situ stresses. 

5. Determine the stress levels within plates and on plate-bounding faults, and therefore 

better understand of the physics of both mantle convection and earthquakes. 

6. Calibrate model rheologies for the lithosphere (both crust and mantle), useful in 

simulations of past and present tectonics at all scales. 

7. Publish a finite element code which, with minor alterations, could be applied to the 

tectonics of Venus, Mars, and the Galilean satellites. 

8. Publish an improved global map of anelastic strain rates and fault slip rates, which will 

yield improved estimates of long-term seismic hazard. 

9. Create a baseline model to estimate ithe results of geodetic campaigns for planning 

purposes, and to highlight the truly novel and unexpected features of real data that 
deserve more intensive investigation. 
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II. Narrative Description of Work 

1992: After receiving our formal notification in April 1992, we began working on the project 
June 1. In 1992, 1 worked one full summer month plus a cumulative total of another month on 
the project. Senior graduate student Xianghong Kong worked full time through the summer 
and half-time during the academic year, for a total of 4.8 months in 1992. There were few 
expenses other than salary, overhead, and my trip to GSFC in October. We had a surplus of 
approximately $15,000 as of the end of 1992, which was partly due to the later-than- 
anticipated start of graduate student support, and partly to the minimal computing expense in 
this initial period of coding and testing. 

The proposal had outlined a year for program conversion, a year for testing and 
debugging, and two years for numerical experiments. We kept to that schedule. In our first 
(partial) year, we designed a finite element for isostatic thin-shell deformation on a sphere, 
derived all of its algebraic and stiffness properties, and embedded it in a new finite element 
code which derives its basic solution strategy (and some critical subroutines) from our earlier 
flat-Earth codes. We also designed and programmed a new fault element to represent faults 
along plate boundaries. We wrote a preliminary version of a spherical graphics program for 
the display of output. We tested this new code for accuracy on individual model plates. We 
made estimates of the computer-time/cost efficiency of the code for whole-earth grids, which 
were reasonable. Finally, we converted an interactive graphical grid-designer program from 
Cartesian to spherical geometry to permit the beginning of serious modeling. 

For reasons of cost efficiency, our models are isostatic, and do not consider the local 
effects of unsupported loads or bending stresses. (To do so would increase computer cost 
by a factor of =16x, and would vastly complicate the algebra. No one has yet attempted this 
with realistic rheologies.) Therefore, we could not employ a standard shell element from the 
engineering literature, and we had to design one. The requirements are: (i) ability to represent 
rigid rotation on a sphere, (ii) ability to represent a spatially uniform strain-rate tensor in the 
limit of small elements, and (iii) continuity of velocity across all element boundaries (the 
intentional discontinuity in fault elements being considered internal). After much trying, we 
concluded that the 6-node isoparametric triangle that we had used in the plane can not be 
adopted to the sphere without violating (iii). Therefore, we designed a 3 -node triangle shell 
element which has two different sets of basis functions to represent (vector) velocity and all. 
other (scalar) variables. Such elements can be shown to converge to the formulas for plane 
triangles in the limit of small size, but can also applied to cover any area smaller than a 
hemisphere. The difficult volume integrals involved in computing the stiffness of such 
elements are performed numerically using 7 Gauss integration points on the surface of the 
sphere, beneath each of which a vertical integral is performed using about 100 points. 

The fault element consistent with these triangles is a simple 4-node line segment. (Its 
dip is pre-assigned.) While such elements cannot elegantly represent the true shape of a 
curving transform fault (like the Altyn Tagh fault), we have programmed them so that the 
direction of slip at each node pair is not necessarily parallel to the element. This means that 
by appropriate assignment of these directions, we can represent strain-free relative rotation of 
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plates along transforms which are arcs of small circles, as well as handle the general case. (In 
elements which are assigned a non-vertical dip, the slip vector has two degrees of freedom, 
and this issue does not arise.) 

The tests of this code in the first year were of two types. First, we verified that the 
new code gives the same result as the old flat-Earth code for small, regional grids in which the 
Earth's curvature is negligible. Second, we made computations on large plates subtending 
about 60®, with velocity boundary conditions. Computing the reaction forces to these 
boundary conditions, we confirmed that torques on these plates are exactly balanced in all 
problems. (If this holds in all cases, it is equivalent to force balance at all points.) 

1993: In the second year (1993), I worked one full summer month plus a cumulative total of 2 
additional months on the project. Senior graduate student Xianghong Kong worked full time 
through the summer and half-time during the academic year, for a total of 8.2 months. 

Besides our salaries and overhead, the principal expenditures were: $2,995 for mainframe 
computing charges, $456 for rental of a color-terminal port on the mainframe, $350 for 
attending the Fall 1992 AGU meeting, $145 for data (NGDC Earth topography CD-ROM), 
$260 for miscellaneous computer supplies (tape backup and ribbons), and $370 for UCLA 
recharges (mail and photographic services). 

The principal accomplishments of the second year were: 

i) Completion of testing and debugging for the main finite element code, SHELLS. As of 
1992, we had only tested the continuum (plate-interior) elements, and not any full models 
including fault elements. These have now been added and debugged. Grids including both 
interior and/or boundary faults can now pass the tests of: 

• essentially zero strain-rate (within numerical precision) under isotropic lithostatic 
loading; 

• essentially zero strain-rate (within numerical precision) under uniform rigid-plate 
rotation about any pole; 

• isotropic horizontal extension (in small plates) under excess topographic load; 

• linear variation of horizontal stresses in horizontal directions under uniform basal 
shear stress; 

• agreement of spherical code with old flat-earth code (PLATES) when the model area 
spans only a few degrees; 

• balance of total torque applied to large plates by boundary conditions, whether 
velocity-type, traction-type, or mixed. 

This tested and verified finite-element program (SHELLS) was presented by Xianghong 
Kong at the Fall 1993 AGU Meeting (T21D-7). 

ii) Programming of an interactive grid-editor program, ORBWEAVER. This program was 
written in QuickBasic for any PC running DOS, with any type of graphics monitor and mouse. 
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It facilitates the creation and editing of finite-element grids (with instant graphical feedback) by 
supporting the following operations: 

• covering a sphere with a uniform grid obtained by subdivision of the icosahedron. 
The subdivision level may be specified, from 0 to 7. Each level has 4 times as 
many spherical triangles as the previous level, and these grids are nested. 

• covering a portion of the sphere (outlined with the mouse) at a subdivision level 
finer than the global level. 

• adding or deleting individual nodes, (triangular continuum) elements, and faults with 
the mouse. 

• dragging nodes to adjust their positions. On fault elements, the compass azimuth at 
the midpoint is continuously displayed for convenience in matching actual fracture 
zone azimuths. 

• testing grid topology for unacceptable gaps, overlaps, or disconnected regions. One 
of the tests used is that the sum of the areas of the spherical triangular elements must 
equal the area within the model perimeter (or total spherical area, if no perimeter). 

• displaying any of the four nodal input data (elevation, heat flow, crustal thickness, 
and mantle lithosphere thickness) in colored contour maps, and permitting 
adjustment of individual values using the mouse and keyboard. 

iii) Creation of program ORBNUMBER, which renumbers the nodes of any grid to minimize 
the bandwidth (which is the maximum difference between the integer numbers in any pair of 
nodes connected by an element or fault). The importance of the bandwidth is that in the main 
finite element calculation, the computer time goes as the square of bandwidth (for a fixed 
number of nodes). This program was written in Fortran?? for the mainframe, and operates by 
assigning all nodes a "connection level" with respect to a starting node. The connection level 
is the minimum number of elements and/or faults intervening between the node and the 
starting node. Beginning with connection level 1, the nodes in the level are numbered in the 
order of their longitude (relative to a pole at the starting node). The resulting sequence of 
node numbers spirals around the globe from pole to antipode and is nearly optimal. 

iv) Creation of a first-generation family of global grids, EARTH2A/B/C.FEG. These grids 
were based on the second subdivision of the icosahedron, modified by hand-editing in 
ORB WE AVER to represent plate-boundary faults. The number of nodes is about 5?2, the 
number of elements is about 750 , and the number of faults is about 185. Frankly, this is quite 
a crude level of resolution; the mean size of the triangular continuum elements is 680,000 
km^, so the mean element side is 11 degrees long. The advantages of a crude grid were: 

• preparation and checking of input data was much faster; and 

• computations were fast and cheap. Specifically, each iteration of any calculation 
with this grid took only 2.5 seconds on the IBM 6090, so a converged model (50 
iterations) took only 2 minutes, and can be computed for $4. This allowed quick 
turnaround and free experimentation with the tactical parameters that affect 
convergence, as well as the physically-meaningful parameters. 
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In drawing the grid, it was necessary to approximate the geometry of spreading rises. 

There are far more transform faults than there are fault elements available. Two solutions 
are possible: (a) use an obliquely-spreading normal fault element to represent many rise 
and transform segments collectively; or (b) keep the orthogonal rise/transform relationship 
by using fewer transforms of exaggerated length. We chose the second approach, because 
only this method preserves the constraint on the direction of the relative velocity vector 
between the two interacting plates. The azimuths of these transforms were carefully set to 
actual azimuths obtained from the NUVEL-1 input data tables of DeMets et al. (1991). 

The topography data for this grid was obtained from the ETOP05 digital dataset. First, 
topography was smoothed with a Gaussian filter of 250 km diameter, to eliminate 
flexurally-supported anisostatic features which our isostatic model cannot simulate. Then, 
topography was sampled at 7 integration points in each element. Finally, a Galerkin 
method was used to find the best-fitting nodal values which (combined with the nodal 
functions) would approximate this surface in a global sense. After this formal procedure, 
we went back and hand-corrected the bathymetry of spreading rises and trenches to their 
actual point values, because the topographic difference between rises and trenches is the 
source of one of the most important plate driving forces. 

Heat flow was sampled from the global dataset of 5 -degree means compiled by Henry 
Pollack of the University of Michigan (personal communication). After the same formal 
procedure, the heat-flow of all rises was set to a uniform 300 mW/m^, to avoid artificial 
strength variations that would otherwise arise from inadequate sampling of hydrothermal 
circulation fields. 

Crustal thickness was calculated from the isostatic assumption, using a lower cutoff of 5 
km. The results have not yet been checked against seismic refraction datasets, but look 
reasonable. 

Total lithosphere thickness was computed by making the base of the mantle lithosphere 
into an isothermal surface. This computation involved a steady-state heat-flow 
assumption, which is not really correct in the oceanic and orogenic regions. This is 
another point for future improvement. 

v) The first real experiments with this global grid were conducted. In the first realistic set of 
models, we imposed the velocity of the subducting slabs as a velocity boundary conditions on 
the footwall side of subduction zone thrust faults. (However, in continental collision belts 
where there may not be any deep slab, we left these footwalls alone.) All velocities were from 
the NUVEL-1 rigid-plate model of DeMets et al. (1991). The base of each plate was 
otherwise free of horizontal traction. 

Different models in this first set differed only in the friction coefficient assigned to plate- 
boundary faults. Low values (e.g., 0.17) were most successful. Velocity results were 
hardly distinguishable from a rigid-plate model. The real information was in the fault- 
slip, intraplate strain-rate, and stress plots. Unfortunately, each contained some probable 
artifacts of the input data, poor resolution, and/or boundary conditions. 
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Many triangular elements lying in the angle between a spreading rise and a transform fault 
had all 3 nodes lying on spreading rises. Therefore, their interiors were at a uniform high 
heat flow equal to the rise value. This made them unreasonably weak, and in many cases 
these elements deformed internally, while the spreading rise either locked, or slipped at 
less than the total plate rate. We corrected this by inserting more and smaller elements 
near the rises and improving the resolution of the input heat-flow model. 

A second serious artifact was the very high value of deviatoric stress in many localities 
near subduction zones. We spent a lot of effort to check that these were not an artifact of 
ill-conditioning of the stiffness matrix, and are satisfied that they are a real part of the 
solution. The problem apparently arises from imposing velocity boundary conditions on 
rather short slab segments, when the global force balance would tend to prefer a different 
plate velocity. This leads to an unreasonable, near-singular concentration of forces. 

These first global models were presented by Peter Bird at the Fall 1993 AGU Meeting 
(T21D-8). 

vi) We also programmed the graphics program used to produce these figures (ORBMAP). 

This program was written in Fortran?? for the mainframe, where we made use of proprietary 
DISSPLA software to drive the output devices. The advantages of DISSPLA were that 
changing map projection is easy, and that the device-independent coding allows us to see the 
plots on a color terminal, then send them to a laser printer (for black-and-white) or to a 
Versatec electrostatic plotter (for color). Program ORBMAP can optionally produce plots of: 

• the finite element grid 

• elevation/bathymetry 

• heat flow 

• crustal thickness 

• temperature of the Moho 

• total lithosphere thickness 

• temperature at the base of the plate 

• any assumed flow pattern below the base of the asthenosphere (which might be 
assumed for a mixed boundary condition on the base of the model) 

• any traction on the base of the plates (from a traction or mixed boundary condition) 

• computed horizontal velocity at the surface 

• velocity changes between iterations (for convergence studies) 

• computed fault slip rates 

• computed intraplate anelastic strain-rates 

• computed (vertically-integrated) deviatoric stress tensors 

All maps produced by ORBMAP are conformal (angle-preserving); global maps are 
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Mercator, and local maps are stereographic. Also, symbols are scaled to take map 
distortion into account, so that any symbol can be compared to the symbol key in the 
margin. 

1994: In the third year (1994), I worked one full summer month (paid) plus about 2 additional 
months (unpaid) on the project. Senior graduate student Xianghong Kong worked full time 
through the summer and half-time during the academic year, for a total of 8.2 months. 

We wrote and submitted a 45 -page manuscript for Journal of Geophysical Research 
detailing our new numerical methods. After a long delay, it was rejected. This was not for 
technical deficiencies, but because it was judged not to be sufficiently interesting for the 
number of pages required! 

In this year, I finished the programming of additional output plot types (velocity- 
discontinuity and most-compressive horizontal principal stress direction). We discovered and 
corrected a minor bug concerning the assumed width of the ductile shear zones underlying 
surface faults in the mantle lithosphere. I did extensive testing to check that no topological 
errors in the grids were impeding plate movement, and coded topological tests into 
ORB WEAVE. We fixed some technical deficiencies in our slab-driven velocity boundary 
conditions, including converting this condition from two components of velocity fixed to only 
one. 


The major technical accomplishment was the switch from subjective to objective 
scoring of models. We selected Robaudo and Harrison (1993) to represent the combined 
results of SLR and VLBI interplate geodesy, and programmed a least-squares adjustment of 
the velocity reference frame between this dataset and the model prior to scoring. Because the 
geodetic data did not cover all plates, we also used tabulated seafloor-spreading rates from the 
NUVEL-1 data tables as an independent test. Horizontal most-compressive stress directions 
from the World Stress Map dataset of Zoback provided a third scalar measure of each model's 
quality (or deficiency). Peter Bird presented these scores at the Fall 1995 AGU meeting. 

In his personal time (50%) and in about half of his NASA-supported time (another 
25%), graduate student Xianghong Kong began a sub-project for his dissertation, which was a 
regional application of our modeling methods to the Eurasian plate. The goals were to provide 
a basis for planning geodetic deployments in the region, to test hypotheses about the eastward 
extrusion (versus shortening) of China, and to determine the rheology of the lithosphere. 

1995: In the fourth year (1995), we were operating under a no-cost extension of time. Kong 
continued to receive half-time support through September. Bird worked a total of about 3 
months on the project without pay. 

Early in the year, we reluctantly revised our previous JGR submission into a 3-page 
technical note (the only option left to us by the editor) and resubmitted this. It was accepted 
promptly in the summer and published in November (see below). AGU agreed to maintain a 
copy of the complete finite-element program in its electronic archives as a supplement to this 
paper. 

Xianghong Kong completed and defended his dissertation on the tectonics of Eurasia. 
He found rapid eastward extrusion (soon to be tested by GPS arrays) and a remarkable 
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weakness (friction < 0.10) of the internal and bounding faults of Eurasia. In his best-fitting 
model, 60% of the strain was by fault slip, and 40% was by block deformation. This result 
tends to retrospectively confirm the appropriateness of our modeling methods which allow for 
both kinds of strain. He and I participated in the Rubey Symposium on the tectonics of 
Eurasia and jointly prepared a chapter for the book (see below). In September 1995, Kong 
took a post-doctoral position and ceased to be paid from this grant. 

The major effort in global modeling this year was a systematic recomputation with the 
final (debugged) version of the finite-element program SHELLS and the final version(s) of the 
EARTH2x finite-element grid. Instead of noisy in-situ heat-flows from the seafloor, we now 
used predicted mean heat-flows based on seafloor age. The results of these 72 experiments 
were fully detailed in our manuscript, "Testing hypotheses on plate-driving forces..." 
submitted last summer, so they will not be repeated here. 

1996: In the past year (1996), the award was no longer active and there were no expenditures. 

I spent about one month in June to port our programs to the new AlX-based computing cluster 
at UCLA, and then prepared the long manuscript describing our 1995 experiments for the 
Journal of Geophysical Research . This preprint was mailed to NASA 22 July 1996. Also, the 
book chapter resulting from Xianghong Kong's dissertation on Eurasia was published (see 
below). 

III. Successes and Failures 

1. Rewrite our existing finite-element code to permit it to be used in global simulations’. 

Done. 

2. Create a forward model of the dynamics of the global lithosphere, permitting 

predictions of how alternative patterns of mantle convection and alternative plate 
rheologies would be expressed as observable geodetic displacements and in situ 
stresses: Done. 

3. Apply this model in a set of systematic simulation experiments to determine these 

quantities: About 200 experiments computed. "Active” mantle convection is 
clearly more successful than "passive" convection. The weakness of plate 
boundary faults found in previous studies is confirmed. However, we have not 
yet found a mantle convection model which is fully satisfactory. The mantle- 
convection community has been challenged to provide one. 

4. First, experiment with parameterized tractions on the bottoms and edges of plates to 

match geodetic plate velocities. Next, optimize tractions on each plate by an iterative 
method: Not attempted. We judged that the results would depend on the 
tractions in such a strongly nonlinear way that it would be prohibitively 
expensive to fully explore the 36-dimensional parameter space. Finally, vary 
plate rheology to match geodetic intraplate strain rates and in situ stresses: Partially 
completed. We optimized the strength of plate-boundary faults. 

5. Determine the stress levels within plates and on plate-bounding faults, and therefore 

better understand of the physics of both mantle convection and earthquakes: We 
confirm the results of previous regional models (CA, AK): Deviatoric stress 
levels are low; topographic stress is a large fraction of the total. Earthquake 
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stress drops are also a large fraction of initial shear stress on these faults. 

6. Calibrate model rheologies for the lithosphere (both crust and mantle), useful in 

simulations of past and present tectonics at all scales*. Not attempted. 

7. Publish a finite element code which, with minor alterations, could be applied to the 

tectonics of Venus, Mars, and the Galilean satellites: Done; code is available by 
Internet or from AGU. 

8. Publish an improved global map of anelastic strain rates and fault slip rates, which will 

yield improved estimates of long-term seismic hazard: Not attempted. The large 
prediction errors remaining in even the best models (13 mm/a, 33°) suggest that 
it would be premature to use them as a basis for seismic hazard estimation. 

9. Create a baseline model to estimate the results of geodetic campaigns for planning 

purposes, and to highlight the truly novel and unexpected features of real data that 
deserve more intensive investigation: Our global model is not yet of this quality. 
However, the regional model of Eurasia serves this purpose. 

IV. Program and Data Products 

Kong, Xianghong, and Peter Bird. "A spherical finite element technique for modeling tectonic 
deformation and fault motion in large regions", EOS Trans. AGU, 74 (43), Fall Meeting 
Suppl., 566 (1993). 

Bird, Peter, and Xianghong Kong. "First global lithosphere models with realistic rheology", 
EOS Trans. AGU , 74 (43), Fall Meeting Suppl., 566 ^993). 

Bird, Peter. "Finite element models of global tectonics require driving forces from mantle 
convection", Eos Trans. AGU , Fall Meeting Supplement, 177 (1994). 

X. Kong and P. Bird: "SHELLS: A thin-shell program for modeling neotectonics of regional 
or global lithosphere with faults," L Geophys. Res., lOQ , no. Bll, 22,129-22,131 (1995). 

X. Kong and P. Bird: "Neotectonics of Asia: Thin-shell finite-element models with faults," in: 
An Yin and T. M. Harrison (ed.s). The Tectonic Evolution of Asia , Cambridge University 
Press, p. 18-34 17 pages (1996). 

Bird, Peter. "Testing hypotheses on plate-driving mechanisms with global models of the 
faulted lithosphere," 36-page manuscript submitted to L Geophys. Res, and to NASA, 7/96. 

All essential programs and datasets produced in this project are available to the public by 
anonymous FTP at: pong.igpp.ucla.edu/pub/pbird/neotec/shells/. 





